*** Set working directory ***
cd "C:\..."

*********************** Load Data **********************************************
clear
graph drop _all
cap drop all

import excel "Data.xlsx", sheet("data") firstrow

************************** Setup ***********************************************
* Choose impulse response horizon
local hmax = 36

* Time variable
gen monthly_date = mofd(date)
format monthly_date %tm
drop date
rename monthly_date date
sort date
tsset date
egen time = group(date)
drop if date == .

* Standardize MP Shock series - unit standard deviation
egen mps_std = sd(chg_sveny02)
gen mps = chg_sveny02/mps_std

sort date
tsset date

* Normalize GSCPI for 1998-2019 sample
rename gscpi gscpi_unscl
egen gscpi_mean = mean(gscpi_unscl)
gen gscpi_dm = gscpi_unscl - gscpi_mean
egen gscpi_std = sd(gscpi_dm)
gen gscpi = gscpi_dm/gscpi_std

* Split mps into pos and neg
gen pos = 0 
gen neg = 0 
replace pos = mps if mps > 0
replace neg = mps if mps < 0

* Interact GSCPI and asymmetric mp shocks
gen gscpi_pos = L.gscpi*pos
gen gscpi_pos1 = L.gscpi*L.pos
gen gscpi_pos2 = L.gscpi*L2.pos
gen gscpi_neg = L.gscpi*neg
gen gscpi_neg1 = L.gscpi*L.neg
gen gscpi_neg2 = L.gscpi*L2.neg

* Put variables in log form
gen cpi_inf = ln(cpi)
gen rgdp_growth = ln(ip)
rename unrt unrt_chg
rename retail retail2
gen retail = ln(retail2)
rename vix vix2

sort date
tsset date

* Interact RHS variables
gen ip_mps = rgdp_growth*mps
gen cpi_mps = cpi_inf*mps
gen unrt_mps = unrt_chg*mps
gen retail_mps = retail*mps
gen ebp_mps = ebp*mps

* LHS variables
forvalues h = 0/`hmax' {
	gen rgdp_growth_`h' = f`h'.rgdp_growth
	gen cpi_inf_`h' = f`h'.cpi_inf
	gen unrt_chg_`h' = f`h'.unrt_chg
	gen retail_`h' = f`h'.retail
}

************************** LP regressions **************************************
/* Run the LPs */
*ssc install lincomest

* Real GDP Growth
eststo clear
cap drop b1 b2 u1a u2a d1a d2a u1 u2 d1 d2 Years1 Zero1
gen Years1 = _n-1 if _n<=`hmax'
gen Zero1 =  0    if _n<=`hmax'
gen b1=0
gen b2=0
gen u1a=0
gen u2a=0
gen d1a=0
gen d2a=0
gen u1=0
gen u2=0
gen d1=0
gen d2=0

forv h = 0/`hmax' {
	* levels
	newey rgdp_growth_`h' l(0/2).pos l(0/2).neg l.gscpi gscpi_pos gscpi_pos1 gscpi_pos2 gscpi_neg gscpi_neg1 gscpi_neg2 l(1/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg l(0/2).retail l(0/2).ebp, lag(`h')
replace b1 = _b[gscpi_pos]                    if _n == `h'+1
replace u1a = _b[gscpi_pos] + 1.645* _se[gscpi_pos]  if _n == `h'+1
replace d1a = _b[gscpi_pos] - 1.645* _se[gscpi_pos]  if _n == `h'+1
replace u1 = _b[gscpi_pos] + 1* _se[gscpi_pos]  if _n == `h'+1
replace d1 = _b[gscpi_pos] - 1* _se[gscpi_pos]  if _n == `h'+1

replace b2 = _b[gscpi_neg]                    if _n == `h'+1
replace u2a = _b[gscpi_neg] + 1.645* _se[gscpi_neg]  if _n == `h'+1
replace d2a = _b[gscpi_neg] - 1.645* _se[gscpi_neg]  if _n == `h'+1
replace u2 = _b[gscpi_neg] + 1* _se[gscpi_neg]  if _n == `h'+1
replace d2 = _b[gscpi_neg] - 1* _se[gscpi_neg]  if _n == `h'+1

eststo
}


* CPI Inflation
cap drop b3 b4 u3a u4a d3a d4a u3 u4 d3 d4 Years2 Zero2
gen Years2 = _n-1 if _n<=`hmax'
gen Zero2 =  0    if _n<=`hmax'
gen b3=0
gen b4=0
gen u3a=0
gen u4a=0
gen d3a=0
gen d4a=0
gen u3=0
gen u4=0
gen d3=0
gen d4=0

forv h = 0/`hmax' {
	* levels
	newey cpi_inf_`h' l(0/2).pos l(0/2).neg l(0/2).gscpi gscpi_pos gscpi_pos1 gscpi_pos2 gscpi_neg gscpi_neg1 gscpi_neg2 l(0/2).rgdp_growth l(1/2).cpi_inf l(0/2).unrt_chg  l(0/2).retail l(0/2).ebp, lag(`h')
replace b3 = _b[gscpi_pos]                    if _n == `h'+1
replace u3a = _b[gscpi_pos] + 1.645* _se[gscpi_pos]  if _n == `h'+1
replace d3a = _b[gscpi_pos] - 1.645* _se[gscpi_pos]  if _n == `h'+1
replace u3 = _b[gscpi_pos] + 1* _se[gscpi_pos]  if _n == `h'+1
replace d3 = _b[gscpi_pos] - 1* _se[gscpi_pos]  if _n == `h'+1

replace b4 = _b[gscpi_neg]                    if _n == `h'+1
replace u4a = _b[gscpi_neg] + 1.645* _se[gscpi_neg]  if _n == `h'+1
replace d4a = _b[gscpi_neg] - 1.645* _se[gscpi_neg]  if _n == `h'+1
replace u4 = _b[gscpi_neg] + 1* _se[gscpi_neg]  if _n == `h'+1
replace d4 = _b[gscpi_neg] - 1* _se[gscpi_neg]  if _n == `h'+1

eststo
}


* Change in Unemployment Rate
cap drop b5 b6 u5a u6a d5a d6a u5 u6 d5 d6 Years3 Zero3
gen Years3 = _n-1 if _n<=`hmax'
gen Zero3 =  0    if _n<=`hmax'
gen b5=0
gen b6=0
gen u5a=0
gen u6a=0
gen d5a=0
gen d6a=0
gen u5=0
gen u6=0
gen d5=0
gen d6=0

forv h = 0/`hmax' {
	* levels
	newey unrt_chg_`h' l(0/2).pos l(0/2).neg l(0/2).gscpi gscpi_pos gscpi_pos1 gscpi_pos2 gscpi_neg gscpi_neg1 gscpi_neg2 l(0/2).rgdp_growth l(0/2).cpi_inf l(1/2).unrt_chg  l(0/2).retail l(0/2).ebp, lag(`h')
replace b5 = _b[gscpi_pos]                    if _n == `h'+1
replace u5a = _b[gscpi_pos] + 1.645* _se[gscpi_pos]  if _n == `h'+1
replace d5a = _b[gscpi_pos] - 1.645* _se[gscpi_pos]  if _n == `h'+1
replace u5 = _b[gscpi_pos] + 1* _se[gscpi_pos]  if _n == `h'+1
replace d5 = _b[gscpi_pos] - 1* _se[gscpi_pos]  if _n == `h'+1

replace b6 = _b[gscpi_neg]                    if _n == `h'+1
replace u6a = _b[gscpi_neg] + 1.645* _se[gscpi_neg]  if _n == `h'+1
replace d6a = _b[gscpi_neg] - 1.645* _se[gscpi_neg]  if _n == `h'+1
replace u6 = _b[gscpi_neg] + 1* _se[gscpi_neg]  if _n == `h'+1
replace d6 = _b[gscpi_neg] - 1* _se[gscpi_neg]  if _n == `h'+1

eststo
}


* Retail Sales
cap drop b7 b8 u7a u8a d7a d8a u7 u8 d7 d8 Years4 Zero4
gen Years4 = _n-1 if _n<=`hmax'
gen Zero4 =  0    if _n<=`hmax'
gen b7=0
gen b8=0
gen u7a=0
gen u8a=0
gen d7a=0
gen d8a=0
gen u7=0
gen u8=0
gen d7=0
gen d8=0

forv h = 0/`hmax' {
	* levels
	newey retail_`h' l(0/2).pos l(0/2).neg l(0/2).gscpi gscpi_pos gscpi_pos1 gscpi_pos2 gscpi_neg gscpi_neg1 gscpi_neg2 l(0/2).rgdp_growth l(0/2).cpi_inf l(0/2).unrt_chg  l(1/2).retail l(0/2).ebp, lag(`h')
replace b7 = _b[gscpi_pos]                    if _n == `h'+1
replace u7a = _b[gscpi_pos] + 1.645* _se[gscpi_pos]  if _n == `h'+1
replace d7a = _b[gscpi_pos] - 1.645* _se[gscpi_pos]  if _n == `h'+1
replace u7 = _b[gscpi_pos] + 1* _se[gscpi_pos]  if _n == `h'+1
replace d7 = _b[gscpi_pos] - 1* _se[gscpi_pos]  if _n == `h'+1

replace b8 = _b[gscpi_neg]                    if _n == `h'+1
replace u8a = _b[gscpi_neg] + 1.645* _se[gscpi_neg]  if _n == `h'+1
replace d8a = _b[gscpi_neg] - 1.645* _se[gscpi_neg]  if _n == `h'+1
replace u8 = _b[gscpi_neg] + 1* _se[gscpi_neg]  if _n == `h'+1
replace d8 = _b[gscpi_neg] - 1* _se[gscpi_neg]  if _n == `h'+1

eststo
}


************************* Plot IRF's *******************************************
* Real GDP
twoway ///
(rarea u1a d1a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u1 d1  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid))  ///
(line b1 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Contractionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_gdp1, replace

twoway ///
(rarea u2a d2a  Years1,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u2 d2  Years1,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b2 Years1, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero1 Years1, lcolor(black)), legend(off) ///
title("Expansionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_gdp2, replace

gr combine fig_gdp1 fig_gdp2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("IP", size(small)) ///

gr rename fig_gdp, replace


* CPI Inflation
twoway ///
(rarea u3a d3a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u3 d3  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b3 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Contractionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_cpi1, replace

twoway ///
(rarea u4a d4a  Years2,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u4 d4  Years2,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b4 Years2, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero2 Years2, lcolor(black)), legend(off) ///
title("Expansionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_cpi2, replace

gr combine fig_cpi1 fig_cpi2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("CPI", size(small)) ///

gr rename fig_cpi, replace


* Unemployment Rate
twoway ///
(rarea u5a d5a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u5 d5  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b5 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Contractionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ur1, replace

twoway ///
(rarea u6a d6a  Years3,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u6 d6  Years3,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b6 Years3, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero3 Years3, lcolor(black)), legend(off) ///
title("Expansionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ur2, replace

gr combine fig_ur1 fig_ur2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("Unemployment", size(small)) ///

gr rename fig_ur, replace


* Retail Sales
twoway ///
(rarea u7a d7a  Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u7 d7  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b7 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Contractionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ebp1, replace

twoway ///
(rarea u8a d8a Years4,  ///
fcolor(gs15) lcolor(gs8) lpattern("-")) ///
(rarea u8 d8  Years4,  ///
fcolor(gs12) lcolor(white) lpattern(solid)) ///
(line b8 Years4, lcolor(blue) ///
lpattern(solid) lwidth(thick)) ///
(line Zero4 Years4, lcolor(black)), legend(off) ///
title("Expansionary MP Shocks", color(black) size(small)) ///
ytitle(" ", size(small)) xtitle(" ", size(small)) ///
xscale(range(0 36)) xlabel(0(6)36) ///
graphregion(color(white)) plotregion(color(white))

gr rename fig_ebp2, replace

gr combine fig_ebp1 fig_ebp2, ///
graphregion(color(white)) plotregion(color(white)) ycommon ///
title("Retail Sales", size(small)) ///

gr rename fig_ebp, replace


* Combine all 4 IRFs
gr combine fig_gdp fig_cpi fig_ur fig_ebp, ///
rows(2) cols(2) graphregion(color(white)) plotregion(color(white)) ///


********************************************************************************
